Again, we can create objects,here start from grid topology definition
topology <- GridTopology(cellcentre.offset = c(1,1,2), cellsize=c(1,1,1), cells.dim = c(3,4,6)) #cellcentre.offset numeric; vector smin coordinates of cell centre in each dim
sp_grid <- SpatialGrid(topology)
class(sp_grid)[1] “SpatialGrid” attr(,“package”) [1] “sp”
Formal class ‘SpatialGrid’ [package “sp”] with 3 slots ..@ grid :Formal class ‘GridTopology’ [package “sp”] with 3 slots .. .. ..@ cellcentre.offset: num [1:3] 1 1 2 .. .. ..@ cellsize : num [1:3] 1 1 1 .. .. ..@ cells.dim : int [1:3] 3 4 6 ..@ bbox : num [1:3, 1:2] 0,5 0,5 1,5 3,5 4,5 7,5 .. ..- attr(*, “dimnames”)=List of 2 .. .. ..$ : NULL .. .. ..$ : chr [1:2] “min” “max” ..@ proj4string:Formal class ‘CRS’ [package “sp”] with 1 slot .. .. ..@ projargs: chr NA
Start to build objects from grid topology, coordinates or read gridded data
x y
1 1 1
2 2 1
3 3 1
4 1 2
5 2 2
6 3 2
7 1 3
8 2 3
9 3 3
10 1 4
11 2 4
12 3 4
[1] "SpatialPoints"
attr(,"package")
[1] "sp"
#from SpatialPoints SpatialPixels
pix_grid <- SpatialPixels(sp_points_grid)
#or this way object gridded (SpatialPoints into SpatialPixels)
class(sp_points_grid)[1] "SpatialPoints"
attr(,"package")
[1] "sp"
[1] FALSE
fullgrid(pix_grid) <- T #set for TRUE and create SpatialGrid
#add atributes (data.frame), to spatial components
atr_grid <- data.frame(matrix(1:12)) #data.frame
atr_grid[1:3,]#first 3 rows of data.frame[1] 1 2 3
[1] "data.frame"
[1] "SpatialGridDataFrame"
attr(,"package")
[1] "sp"
Visaualisation with lattice graphics
Common way of importing gridded data; calculate summary statistics…
dem_1k.asc has GDAL driver FITS and has 21 rows and 21 columns
Min. 1st Qu. Median Mean 3rd Qu. Max.
0,000 6,419 157,996 301,073 511,735 1255,097
As example plot histogram of elevation values
ESRI asc, SAGA sgrd…
Extremely popular for analysis of big rasters
Data does not need to be in memory
Easy shift to ‘sp’ classes
Extract wanted spatial extent or aggregate by polygon, select on points…
Values of Digital elevation model for each county
dem_1k.asc has GDAL driver FITS and has 21 rows and 21 columns
dem_r <- raster(dem)
#extracted values for each 1 km cell in polygons
#object of class list, length 21, as many counties
dem_county <- extract(dem_r, HRV_adm1_ETRS)
dem_county[1:20][[1]] [1] 722,6663 261,4872 370,5971 773,7128 121,7043 307,7529 601,9179 259,0665
[[2]] [1] 157,0298 122,3399 136,0429 153,9509 157,9960 149,9645 356,1018
[[3]] [1] 125,55239 143,63921 144,06564 97,78848
[[4]] [1] 98,8438 387,1294
[[5]] [1] 196,9077 168,1451 122,0084
[[6]] [1] 185,08521 377,42798 163,93018 281,47314 64,61814 173,73517 24,62883
[[7]] [1] 145,2687 143,5296 336,7780 201,0698 206,9406 757,8558 546,8431 [8] 352,1976 266,6391 691,3894
[[8]] [1] 198,0352 150,9576 123,2260 138,8395 153,0450
[[9]] [1] 285,3583 284,5465 202,5881 237,4919
[[10]] [1] 524,9102 652,7065 429,8729 710,6260 729,8898 779,5086 760,5005 [8] 746,7065 675,1872 809,3881 826,6930 702,2512 725,9460
[[11]] NULL
[[12]] [1] 92,74722 90,60442 84,34888 80,32372 142,12306 95,99942 78,94941 [8] 115,79662
[[13]] [1] 176,4084 382,2073 293,2750 217,3138 196,9911
[[14]] [1] 964,98822 497,56903 457,62909 807,84265 762,32831 171,38272 386,63550 [8] 97,36595 48,06548 26,26266
[[15]] [1] 95,6325 126,6646 141,5341 131,3787 105,1319 111,2801 188,6825 [8] 269,6385 188,2933 125,9783 133,8734 272,3000
[[16]] [1] 714,4169 371,5448 444,0142 626,7087 280,6295 547,9183 355,3295
[[17]]
228,8269
[[18]] [1] 114,6883 191,4677 153,4848 103,0313 412,7148 267,1815
[[19]] [1] 81,57622 75,48589
[[20]] [1] 831,36621 800,62671 16,76544 168,52444 397,27487 534,42041 62,76883 [8] 200,08041 56,51281
Ask for maximal value of elevation by county
Warning in FUN(X[[i]], ...): no non-missing arguments to max; returning -
Inf
[[1]] [1] 773,7128
[[2]] [1] 356,1018
[[3]] [1] 144,0656
[[4]] [1] 387,1294
First we need to prepare data
#import geo data into R
library(rgdal)
#read via GDAL/OGR
#http://thematicmapping.org/downloads/world_borders.php
W <- readOGR(getwd(),"TM_WORLD_BORDERS-0.3", verbose = F, stringsAsFactors=FALSE)
#or
#W <- download.file("http://www.mappinghacks.com/data/TM_WORLD_BORDERS-0.3",getwd(), method="auto",mode="wb")CRS arguments: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
FIPS ISO2 ISO3 UN NAME 0 AC AG ATG 28 Antigua and Barbuda 1 AG DZ DZA 12 Algeria 2 AJ AZ AZE 31 Azerbaijan
Zoom, additional improvements
#create gridlines
grat <- gridlines(W, easts=seq(-40,20,by=5), norths=seq(10,80,by=5), ndiscr = 20)
#zoom to specified spatial extent on plot
plot(W, xlim=c(-30,30), ylim=c(30,58),col='olivedrab3',bg='lightblue')
lines(grat, lty=2)#add gridlines
title(main="My map with zoom", sub="Gridlines added") #add titleload("df.RData") #load data.frame saved in R format
#create spatial objects from data.frame
coordinates(df) <- ~ lng + lat
#define CRS to the points
#World Geodetic System - WGS
proj4string(df) <- CRS("+proj=longlat +datum=WGS84") #EPSG:4326
#colors(2)
col_2 <- colors()[351]
x <-Polygon(cbind(x=c(13.9, 20.2, 20.2, 13.9), y=c(42.5, 42.5,46.5, 46.5)))
class(x)[1] “Polygon” attr(,“package”) [1] “sp”
Formal class ‘Polygon’ [package “sp”] with 5 slots ..@ labpt : num [1:2] 17 44,5 ..@ area : num 25,2 ..@ hole : logi TRUE ..@ ringDir: int -1 ..@ coords : num [1:5, 1:2] 13,9 20,2 20,2 13,9 13,9 42,5 42,5 46,5 46,5 42,5
plot(W, col=col_2, xlim=c(13.9, 19.6), ylim=c(42.5, 46.5), main="Hello Croatia", lwd=2, axes=T, xlab="Longitude", ylab="Latitude", col='lightgray',bg='lightblue')
#add points to existing graph
points(df, pch=21, cex=1.75, col="red4")
axis(1, col.axis = "dark gray")
axis(2, col.axis = "dark gray")
box(lty = 'solid', col = 'light gray')
text(x=14.7, y=46, labels="AUSTRIA", cex=0.75)
text(x=17.8, y=44.5, labels="BOSNIA AND HERZEGOVINA", cex=0.75)
text(x=18, y=46.4, labels="HUNGARY", cex=0.75)
text(x=19.2, y=42.7, labels="MONTENEGRO", cex=0.75)